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INTRODUCTION 


This  is  the  final  report  of  a  multi-year  research  project  for  modeling 
vortical  flows  with  adaptive  Navier-Stokes  algorithms.  The  objectives  of 
the  research  were  to: 

Develop  a  three-dimensional  adaptive  Navier-Stokes  computational 
model, 

Investigate  adaptation  criteria  and  turbulence  models  suitable  for 
vortical  flows,  and 

Compute  delta  wing  and/or  canard-wing  flows,  analyze  the  results,  and 
compare  with  data. 

Considerable  progress  was  made  in  understanding  the  physics  of  vortical 
flows  through  three-dimensional  calculations,  interactive  flow 
visualization,  and  theoretical  modeling.  Of  particular  interest  was  vortex 
bursting,  both  steady  and  non-steady.  An  adaptive  Euler  equation  solver 
was  developed  and  utilized  to  calculate  flows  for  simple  delta  wing 
geometries  as  well  as  an  F117A  configuration.  The  objective  of  developing 
adaptive  Navier-Stokes  capability  was  not  achieved.  This  was  partly  due  to 
the  fact  that  it  took  longer  than  we  anticipated  to  develop  the  Euler  part  of 
the  solver,  and  partly  to  the  fact  that  vortex  bursting  is  sensitive  to  viscous 
effects  for  sharp  leading  edge  geometries.  During  the  grant  period  the 
Principal  Investigator  became  aware  of  the  interesting  and  important 
aerodynamic  problem  of  normal  force  hysteresis  for  rapidly  pitching  wings 
and  aircraft.  The  hysteresis  is  caused  by  a  lag  in  the  vortex  bursting 
position.  Computational  studies  were  undertaken  to  understand  this 
problem,  even  though  unsteady  flows  were  not  in  the  original  list  of 
objectives. 

A  number  of  publications  emerged  from  the  research,  and  all  are  included 
in  the  list  of  references.  In  the  following  pages,  the  most  important 
findings  from  this  work  are  summarized.  The  reader  is  referred  to  the 
individual  publications  for  further  details. 
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RESEARCH  RESULTS 


Adaptive  Euler  Equation  Solver 

A  considerable  portion  of  the  three  year  grant  was  focused  on  the 
development  of  a  code  for  an  adaptive  unstructured  (tetrahedra)  grid 
solver  for  the  Euler  equations.  The  underlying  algorithms  were  developed 
and  implemented  from  first  principles.  A  trilinear  Galerkin  Finite  Element 
Method  was  used  for  spatial  discretization  and  the  Jameson  four  stage 
Runge-Kutta  type  scheme  was  used  for  temporal  integration.  A  blended 
second-fourth  difference  artificial  viscosity  based  upon  the  Holmes- 
Connell-Saxer  approach  was  adopted.  Boundary  conditions  for  far  field, 
solid  body,  symmetry  plane  and  sharp  edges  were  implemented.  This 
explicit  method  was  run  using  local  time  steps  for  steady  state  problems. 
For  unsteady  problems,  a  novel  acceleration  factor  was  introduced  which 
provided  an  order  of  magnitude  reduction  in  computing  time  for  a  three- 
dimensional  problem  (Reference  11).  A  global  time  step  was  chosen  and 
used  for  all  cells  where  the  minimum  time  step  exceeded  the  global  value. 
For  the  other  cells,  the  local  time  step  was  used.  Although  this  approach 
sacrifices  locally  accuracy  in  some  cells,  a  priori  and  computational 
experiments  established  acceptable  levels  of  global  accuracy.  The  method 
is  very  robust  and  extremely  easy  to  implement.  It  is  ideal  for  parallel 
applications  avoiding  the  difficulties  of  block  factored  implicit  methods.  A 
full  account  of  the  solver  is  given  in  Reference  81  and  an  abbreviated 
account  is  given  in  Reference  12. 

Adaptation  parameters  for  vortical  flows 

In  results  reported  in  Reference  2,  we  experimented  with  different  criteria 
for  adaptive  calculations  of  vortical  flows.  Most  adaptive  calculations  we 
are  aware  of  have  been  designed  to  detect  and  refine  shock  waves.  One 
cannot  expect  the  same  feature  finders  used  for  shock  waves  will  work 
with  vortices.  We  explored  two  parameters  for  adaptation.  One  was  the 
total  pressure  loss  which  is  well  known  to  occur  in  vortical  flow.  The  other 
was  normalized  helicity,  the  dot  product  of  vorticity  and  velocity 
normalized  by  the  absolute  values  of  the  two  quantities.  This  is  just  the 
cosine  of  the  angle  between  the  velocity  and  vorticity  vector  which  should 
be  unity  in  the  core  of  a  vortex.  A  full  discussion  of  the  adaptation 
methods  is  given  in  Landsberg's  thesis.  Reference  2.  The  key  results  are 
contained  in  a  review  paper  written  by  Landsberg  and  Murman  (Reference 


1  A  listing  of  the  code  is  available  upon  request. 
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7)  which  discusses  various  methods,  including  adaptation,  for  controlling 
the  numerical  diffusion  of  vortices. 

In  the  extension  to  unsteady  flows  reported  in  References  8  and  12, 
entropy  was  used  instead  of  total  pressure.  Entropy  is  a  state  variable 
independent  of  the  reference  frame.  Our  conclusion  was  that  entropy  is 
easy  to  implement  and  a  suitable  parameter  for  finding  vortices,  at  least 
for  the  Mach  number  range  considered  in  this  study. 

Vortex  Breakdown  Calculations 

The  adaptive  Euler  Solver  summarized  above  was  used  to  compute  the 
flow  about  a  delta  wing  with  75  degree  leading  edge  sweep,  a  Mach 
number  of  0.3,  and  for  angles  of  attack  from  0  to  52  degrees.  The  base 
non-adapted  grid  had  about  15,000  nodes  while  the  adapted  grid  had 
about  80,000  nodes.  Normal  force  results  agree  well  with  experimental 
data  through  the  entire  range  of  angles  of  attack  when  adaptation  is  used. 
Without  adaptation,  the  grid  resolution  is  insufficient  to  capture  vortex 
breakdown.  Results  for  pitching  moment  and  vortex  breakdown  position 
are  in  reasonable  agreement  with  the  data.  Computations  with  breakdown 
present  do  not  reach  a  steady  state  solution,  most  probably  due  to  the 
unsteady  nature  of  the  burst  vortex.  Small  oscillations  in  forces  and 
bursting  position  persist.  However,  the  solutions  are  steady  in  the  large. 
Comparisons  at  several  angles  of  attack  are  made  with  published  results 
using  a  much  finer  non  adaptive  structured  grid  and  agreement  is 
satisfactory.  Detailed  analysis  of  the  flow  field  reveal  interesting  features 
of  the  breakdown.  At  32  degrees  angle  of  attack,  a  spiral  breakdown  is 
predicted.  At  42  degrees  angle  of  attack,  a  bubble  breakdown  followed  by 
reformation  of  the  vortex  and  a  spiral  breakdown  is  predicted.  Detailed 
results  are  given  in  References  8  and  12. 

Calculations  for  the  F117A 

The  adaptive  Euler  solver  was  applied  to  the  F117A  (Reference  13) 
configuration  to  test  its  ability  to  handle  arbitrary  geometry.  A  plastic 
model  was  digitized  for  the  geometrical  definition.  Results  were  obtained 
at  15,  20,  and  25  degrees  angle  of  attack.  Vortex  bursting  was  absent  in 
the  first  case.  For  20  degrees  angle  of  attack,  bursting  was  about  at  the 
trailing  edge.  Bursting  was  well  over  the  wing  for  the  third  case. 

Subsequent  wind  tunnel  studies  have  shown  general  agreement  with  these 
findings. 
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During  the  grant  period,  we  were  made  aware  of  the  interesting 
experimental  results  which  show  significant  normal  force  hysteresis  for 
pitching  delta  wings  and  for  a  Russian  aircraft.  Dynamic  force  coefficients 
can  exceed  static  values  by  30-50%  on  the  upstroke  portion  of  the  pitch, 
and  can  undershoot  static  values  by  lesser,  but  still  significant  values  on 
the  down  stroke  portion.  Experimental  evidence  shows  this  is  due  to  a  lag 
in  the  vortex  bursting  position,  thereby  maintaining  the  augmented  normal 
forces  to  a  higher  angle  of  attack  on  the  upstroke.  A  lag  in  the  opposite 
direction  on  the  down  stroke  causes  the  undershoot.  A  computational 
experiment  was  done  to  determine  if  this  effect  could  be  predicted  by  the 
adaptive  Euler  solver.  Preliminary  results  were  obtained  (Reference  12) 
showing  the  general  effect.  Both  time  and  resources  ran  out  before  more 
detailed  calculations  could  be  done. 

Visualization  of  vortical  flows  on  unstructured  grids 

Visualization  of  the  complex  vortical  flows  proved  to  be  an  interesting  and 
important  challenge  in  order  to  understand  the  results  (References  1,  3,  4, 
5,  6).  Considerable  attention  was  paid  to  developing  novel  visualization 
strategies  and  evaluating  their  usefulness.  The  foundation  for  the 
visualization  work  is  the  VISUAL  3  software  developed  by  Giles  and 
Haimes  in  the  MIT  Computational  Fluid  Dynamics  Laboratory.  This 
software  treats  generic  three-dimensional  data  sets  on  unstructured  grids, 
and  it  runs  interactively  on  high  end  graphics  workstations.  To  VISUAL  3 
we  added  a  number  of  capabilities  which  provide  for  a  hierarchy  of 
approaches  for  analyzing  complex  flow  fields.  The  strategy  is  to  start  with 
Identification  methods  which  locate  the  major  regions  of  interest  in  the 
flow  field.  Several  new  identification  methods  were  developed  including  X- 
rays,  vector  clouds,  and  a  shock  finder.  When  the  important  regions  are 
found,  the  next  step  is  to  use  Scanning  methods  to  zero  in  on  detailed 
features.  Cutting  planes  and  iso  surfaces  are  examples  of  these.  The  final 
step  is  to  extract  quantitative  data  with  Probing  techniques.  A  variety  of 
interactive  probes  were  developed,  and  significantly  improved  particle 
path  methods  were  introduced.  The  latter  permit  viewing  not  only  the 
path  tangent  to  a  vector  field,  but  also  its  divergence  and  curl  to  be  seen. 
Methods  for  determining  grid  quality  were  also  investigated  (Ref  9). 

Modeling  of  Vortex  Breakdown 


In  order  to  better  understand  the  ability  of  inviscid  computational 
methods  to  capture  vortex  breakdown,  studies  were  undertaken  to 
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compute  and  analyze  an  isolated  axisymmetric  vortex  subjected  to  a 
controlled  adverse  pressure  gradient.  An  incompressible  Navier-Stokes 
solver  was  written  and  applied  to  a  vortex  in  a  tube  with  a  converging- 
diverging  nozzle.  The  breakdown  location  can  be  controlled  by  the  inlet 
swirl,  the  Reynolds  number,  and  the  divergence  angle  of  the  nozzle.  Of 
these  parameters,  the  Reynolds  number  effects  are  insignificant  above 
relatively  low  values. 

Theoretical  studies  have  also  been  done  from  the  point  of  view  of  vorticity 
dynamics.  It  was  found  the  important  effect  leading  to  vortex  breakdown 
is  the  tilting  of  the  axial  vorticity  vector  into  the  azimuthal  direction  as  the 
streamtubes  expand  radially  due  to  the  adverse  pressure  gradient.  This  is 
an  inviscid  effect  which  is  in  agreement  with  the  observation  that 
Reynolds  number  plays  an  insignificant  role.  These  results  will  be 
presented  at  an  upcoming  AIAA  meeting  (Ref  10). 

Turbulence  modeling 

For  Reynolds  numbers  applicable  to  flight,  leading  edge  vortices  have 
turbulent  cores.  Yet  little  work  has  been  done  to  develop  turbulence 
models  for  use  in  CFD  studies  of  vortices.  Most  investigators  who  solve  the 
Reynolds  Averaged  Navier-Stokes  equations  use  turbulence  models 
developed  for  wall  bounded  flows.  There  is  no  reason  to  believe  these  are 
suitable  for  vortical  flows.  Preliminary  work  was  done  to  develop  a 
turbulence  model  for  vortex  flows.  The  results  are  given  in  the  Appendix. 
However,  when  the  goal  of  developing  a  Navier-Stokes  solver  was 
abandoned,  this  task  became  of  lesser  importance  and  it  was  abandoned. 

Parallel  processing  study 

We  expended  some  effort  on  mapping  the  above  unstructured  Euler 
program  onto  the  Kendall  Square  Research  parallel  computer  prior  to  the 
date  it  was  available.  Our  studies  showed  that  the  unstructured  grid 
program  would  achieve  about  6  gigaflops  on  the  full  IK  node  processor. 
Although  wu  expected  to  implement  the  results  on  the  KSR  machine,  delays 
in  its  avaiUMlity  prevented  us  completing  this  task. 
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APPENDIX 


TURBULENCE  MODEL  FOR  A  ROTATING  FLOW 
M.T.  Landahl 


1.  Introduction 

The  modelling  of  the  turbulent  stresses  in  a  flow  with  strong  rotation  is  an  important  problem  in  the 
calculation  of  flow  in  rotating  machinery  and  in  strong  vortices.  How  to  include  the  effect  of 

rotation  in  turbulence  models  like  the  k -e-model  has  been  discussed  by  many,  see  Speciale  !991). 
A  recent  treatment  of  the  problem  by  Ekander  and  Johansson  (1989)  has  led  to  a  simple  algebraic 
Reynolds  stress  model  (ARSM)  which  appears  to  include  the  effect  of  rotation  in  a  consistent 
manner. 

2.  Model  fundamentals 

We  will  consider  an  incompressible  turbulent  flow  in  a  Cartesian  frame  of  reference  rotating 

steadily  around  the  z-axis  by  the  angular  velocity  w.  We  subdivide  the  velocity  components  and 
pressure  in  mean  and  fluctuating  components  by  setting 


Ui(xj.t)  -  Uj  +  uj(xj,t),  P  =  P  +  p(xj.t),  <u{>=0,  <p>=0  ( 1  ) 

The  momentum  equations  for  the  mean  quantities  (assumed  steady)  read 


where  the  Reynolds  stresses  t  jj  are  defined  by 

t  jj-  -  r<ujUj>  (3) 


and  P*  is  the  reduced  pressure  including  the  centrifugal  effect.  The  fluctuating  components  satisfy 


and  where 


D  _  d  .  fj.  d 

a=s+ui5i; 

tjj  =  -r(ujuj  -  <ujUj>) 


(4) 

(5) 


(6) 

(7) 
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For  the  two-dimensional  case  with  rotation  araound  the  z-axis  the  equations  for  the  fluctuating 
components  read,  with  ui=  u,  U2  =  v,  U3  =  w 


Du  au  ^  au  i-V  a,mau  t  . , 

Br  +  U5T  +  V3y  ‘ 2wv-  ■rt'ST 

(8) 

i? + “  Jr + '$* 2wu-  *ij>  1 

(9) 

IT* 

(10) 

The  equations  for  the  Reynolds  stresses  of  interest,  1 1 1 «  -  r<u^>  1 12=  -r<uv>  and  1 22“  *  r<v2> 
are  obtained  by  multiplying  the  equation  for  u  by  u,  then  by  v,  and  that  for  v  by  u,  then  by  v  and 
averaging  and  adding.  After  using  the  equation  of  continuity  (5)  (me  obtains  transport  equations  for 
the  Reynolds  stresses  which  may  be  written,  following  Ekander  and  Johansson  (1989)  as  follows: 

-jyf1 '  D,i  ~  p|i +  p,i +  e,i 

(11) 

where  is  the  "true"  time  rate  of  change  defined  by 

D’t  “  Dt  +E‘J 

(12) 

with  Eij  being  the  redistribution  of  the  Reynolds  stresses  due  to  the  rotation,  given  in  the  two- 
dimensional  case  by 

Ejj  =  -2w  <uv>,  E12  =w  (<u2>  -<v2>),  E22  =2w  <uv> 

(13) 

Djj  is  the  diffusive  rate  of  change  defined  by 

D  8<ukUjUj>  ld<Ujp>  l8<Ujp> 

^  3xk  r  3xj  r  3xj 

(14) 

Pjj  the  production 

m  -  dUi  -  dUi 

w0ijkt  ik  -  weijktjk 

(15) 

Pjj  the  pressure-strain  correlation 

o  1  /  dul  1  ,  dm 

P|i  T* 

(16) 

and  0jj  the  viscous  dissipation, 
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(17) 


_  „  3ui  .9uj  dui 

e9  '  ‘  n<3^  <'-1+ 


_!<,>  n<^l 
3xk  '3xk'dxi'  3xk 


It  should  be  noted  that  one  half  of  the  contribution  resulting  from  the  Coriolis  term  has  been  kept 
on  the  left-hand  side  of  the  evolution  equation  (8)  on  the  argument  that  it  should  be  consideresd 
part  of  the  substantial  rate  of  change.  This  results  from  the  requirement  that  the  evolution  equations 
should  be  invariant  under  arbitrary  Galilean  transformations.  A  problem  with  earlier  treatments  of 
the  problem  has  been  that  all  of  this  term  has  been  incorporated  incorrecdy  into  the  production 
term. 


By  taking  the  trace  of  the  Reynolds  stresses  one  obtains  the  usual  equation  for  the  kinetic 
energy:  with  k  =  (l/2)<ujUj> 


§r(l/2r)lxi-  d/2)(Pii+  Pii-Djj)-  e  (18) 

where  e  is  the  dissipation  rate 


e  =  2n<^i+|^L)2>  (19) 

oX j  dXj 

For  the  two-dimensional  case  considered  one  thus  finds  (note  that  Ejj-O) 

a=.<uv>^-e  (20) 

For  the  individual  production  terms  P,j  we  have  in  this  case 


Pjl  =  -  2<uv>^  -2w<uv>,  P12  =  ■<v2>^  +  w(<u2>  -  <v2>),  P22  =  2w<uv>(21) 

In  the  case  of  a  curvlinear  coordinate  system  the  derivatives  of  the  base  vectors  will  enter  in  the 
redistribution  terms  Ejj .  For  a  circular  channel  with  rotation  about  the  origin  one  finds  with  uj  =  ur 

,  U2  =  uq  ,  and  t  j  j*  -  r<ur  ur  >,  1 12“  -  r<ur  uq  >,  1 22*  *  rcuquq  > 


Pll  =  2<uruq>^ ,  P12  =  -  2<UfUr>^42<uqUq:>^,  P22  =  -  2<uruq>^q  (22) 
En  =  2<ufuq>^  ,  E12  =(<uquq>-<urur>)  ^9  e22  =  -  2<uruq>^a  (23) 


3.Cons traction  of  an  algebraic  stress  model 


The  fundamental  idea  is  to  construct  algebraic  reladonsships  between  the  individual  components 

t  y  and  the  kinetic  energy  k  and  the  dissipation  rate  6.  To  this  purpose  an  approximation  suggested 
by  Rodi  (1976)  was  utilized  in  which  he  suggested  that  the  difference  between  the  advective  and 

diffusive  terms  of  the  individual  components  t  jj ,  divided  by  the  corresponding  difference  in  the 
equation  for  k,  should  be  set  equal  to  the  ratio  between  the  respective  Reynolds  stress  component 
and  k.This  then  leads  to  the  relation 


D’tii 

D't 


Dk 

Dt 


Dm 


(24) 


In  this  manner  the  problem  of  detemriniing  the  Reynolds  stresses  is  simplified  to  that  of  solving  the 
evolution  equations  for  k  and  6.  For  these  one  needs  modelling  of  all  the  terms,  except  the 
production  and  the  redistribution  terms. 


The  'slow  pan'  of  die  pressure  strain  term  (i.e.,  the  pan  involving  nonlinear  fluctuation  terms  in 
the  pressure  equation)  is  modeled,  following  Rotta  (1951)  as  proportional  to  the  anisotropy  of  the 
Reynolds  stress  tensor : 

pijs  -  ■  fkdij)  (25) 

where  Cg  =  1.9.  The  rapid  pan  is  modeled,  following  Reynolds  (1970),  as  being  proportional  to 
the  anisotropy  of  the  production  tensor 

pi|r-  (26) 

where  P  is  the  trace  of  the  production  tensor  and  c,  =  0.6.The  presence  of  die  wall  will  also 
transfer  energy  from  the  transverse  to  the  stream  wise  component  through  reflections.at  the  wall.  In 
the  paper  of  Ekander  and  Johansson  (1989)  this  effect  is  modelled  in  a  simple  manner  by 
subtracting  bom  the  noimaal  component  die  value 

-cwf<Paa-§P)L  (27) 

and  adding  this  to  the  stream  wise  component  In  (27),  *  0.72  and  L  =  min  (Vy,  1.0),  y 

denoting  the  dHhncr.  to  the  wall  and  1$*  (C|k)3^/(©k)  with  Cf«0.24  and  k=0.41  (the  von 
Karman  constant).  Also,  the  Greek  index  symbols  are  used  to  indicate  that  summation  shall  not  be 
applied.  The  diffusion  rate  is  modelled,  following  Daly  and  Harlow  (1970),  as 

D  *  *  0.22  (28) 

For  the  dissipation  it  is  usually  assumed  that  the  smallest  scales,  which  presumably  carry  most  of 
the  viscous  dissipation,  are  isotropic,  hence 
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eij  =  dij  2e/3  (29) 

The  standard  modelled  equation  for  e  is  modified  to  be  consistent  with  the  k -equation  as  follows: 


with  the  following  constants  chosen:  ce  =  0.18,  cie  =  1.44 ,  C2  e  =  1.92. 

In  their  application  to  channel  flows  Ekander  and  Johansson  (1989)  also  had  to  apply  a  wall 
correction  for  the  inner  part  of  the  inner  inertial  layer,  y*  <  30.  They  used  a  form  proposed  by 
Rodi  (1980). 

Their  results  for  the  flow  in  a  nonrotating  channel,  a  rotating  two-dimensional  channel,  and  for  a 
circular  channel  showed  good  both  qualitative  and  quantitative  agreement  with  experiments.  Their 
results  are  reproduced  in  Figs.  1—7. 

4,  Conclusions 

The  model  proposed  by  Ekander  and  Johansson  (1989)  for  turbulence  in  a  rotating  frame  of 
reference  appears  self-consistent  and  gives  results  in  good  agreement  with  experiments.  Their 
model  should  be  applicable  to  die  flow  in  a  strong  vortex. 
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